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Abstract 

This article describes a bridge between POD-based model order reduction techniques 
and the classical Newton/Krylov solvers. This bridge is used to derive an efficient algo- 
rithm to correct, "on-the-fly" , the reduced order modelling of highly nonlinear problems 
undergoing strong topological changes. Damage initiation problems are addressed and 
tackle via a corrected hyperreduction method. It is shown that the relevancy of reduced 
order model can be significantly improved with reasonable additional costs when using 
this algorithm, even when strong topological changes are involved. 

Keywords: Model order reduction (MOR); Time-dependant nonlinear mechanical prob- 
lems; Proper orthogonal decomposition (POD); Newton/Krylov Solver; Projected conju- 
gate gradient; System approximation; Hyperreduction; CH-POD; Damage propagation; 



1 Introduction 

The simulation of failure in complex materials has been one important issue throughout en- 
gineering. Micro and mesoscale models have demonstrated their ability to predict complex 
phenomena such as crack initiation and propagation in various fields: composite materials 
in aeronautics [U [2] , cement-based structures [3] in civil engineering or biostructures [3] in 
medical engineering. However, these models usually describe the material at a very fine 
scale compared to the size of the structure, and the finite element discretization of the un- 
derlying partial differential equations leads to large, potentially highly nonlinear (therefore 
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requiring a fine discretization in time), numerical problems. Various multiscale computational 
strategies have been developed to tackle this important issue, such as enrichment techniques 
[5J E] , homogenization strategies [3 [H] , domain decomposition methods |TU] , or model order 
reduction [Til [12]. The most efficient of these techniques couple several of these advanced tools 
to efficiently perform multiscale simulations [131 E3 IS Ell [TTJ, [18] . 

Among them, Proper Orthogonal Decomposition-based (POD) [11 ). 112 ^ [T9] model order re- 
duction strategies provide extremely valuable tools for an automatic derivation of a multiscale 
computational method (no a priori physical understanding of the different scales involved is 
required). They consist in using a set of potential solutions to the initial problem (snapshot) 
and extracting, by a spectral analysis a few basis vectors spanning a space of small dimension 
in which the solution to the, initially large, numerical problem is well approximated [20 1 . 
These techniques have been proved extremely efficient when speed is more important that 
accuracy, for instance if the goal is to provide an approximation of the response of a com- 
plex structure in real-time (an example of biomedical application is given in [21]), or quasi 
real-time (in the case of early stage design). This type of applications is our main focus. 

In the context of structural problems involving plasticity or damage, two severe drawbacks 
limit the direct application of POD-based model order reduction: 

• The initial snapshot might be too poor to represent accurately the solution of the dam- 
aged structure. This can happen if: 

— The "a priori" mechanical understanding of the structure is too poor, which might 
result in a bad choice of the snapshot simulations. 

— The number of parameters involved is too important. The snapshot needed to 
compute a relevant reduced basis might be very large, and its reduction inefficient. 

— Strong topological changes in the structure occur. Improvements have been pro- 
posed by Ryckelynk's team [22J where the reduced basis is enriched during the 
computation. When required, a Krylov subspace associated to a linearised op- 
erator of the initial problem is generated, and a few additional basis vectors are 
obtained by a spectral analysis of this space. However, it seems that, when deal- 
ing with complex constitutive laws, this adaptive strategy can be computationally 
expensive [T8] . 

• The integration of the constitutive law needs to be done at each integration point, 
regardless of the dimension of the reduced space. These problems have been handled in 
|23j by a technique called hyperreduction, which consists in considering only a few local 
residuals to compute the internal forces. 

On the other hand, various model order reduction strategies have also been used to de- 
rive good initializations and preconditioners for iterative algorithms for complex nonlinear 
problems at the fine scale |23l I55| I5| H5| ITB| . 

We propose here a novel strategy that couples these two approaches. Our vision is that if 
complex changes in the topology of the structure appear (local crack initiation for instance), 
they can only be accurately predicted by solving the fine scale model, at least locally. This 
can be done, for instance, by adding the finite element shape functions to the snapshot, 
or by using relocalization techniques. However, the long-range effects of these topological 
changes do not require the fine description, and can be obtained by solving the full system 
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very crudely. Hence, when the global residual exceeds a given threshold value, we propose to 
perform conjugate gradient iterations on the full system, orthogonally to the snapshot space 
by means of a classical projection. The new vector obtained by the conjugate gradient does not 
belong to the initial approximation space and is added to the reduced basis. It is interpreted 
as an "on-the-fly" enrichment of the reduced model. The proposed algorithm largely differs 
from what had been done in [22J in the sense that only a single vector belonging to the Krylov 
algorithm is used to enrich the reduced basis, for each of the corrections performed within a 
time step. It thus keeps the size of the reduced basis to a minimum. In addition, the solution 
of the current time step is not restarted after enriching the reduced basis. At last, using 
the properties of the projected Krylov during the correction steps automatically ensures the 
required orthogonality between the additional reduced basis vectors and the initial reduced 
basis. 

This strategy can be interpreted as a bridge between "exact" Newton-Krylov strategies, 
and projection-based model order reduction methods (like the POD) for nonlinear problems. 
Indeed, setting the residual threshold to a low value leads to the former strategy, while setting 
it to a high value leads to the latter. An intermediate value yields an adaptive model order 
reduction method with control of the global residual. We will show that this method permits 
to obtain a correction of the reduced model at cheap price, when the current reduced basis is 
no more sufficient to represent the solution accurately. 

The paper is organized as follows. In section 1, the nonlinear system of equations resulting 
from the discretization of a continuum mechanics problem is introduced, and the model order 
reduction and associated time solution procedure are described. In section 2, we perform the 
reduction of a damage model, and show that the basic application of the proper orthogonal 
decomposition is inadequate. The corrective reduced order modelling is introduced in sec- 
tion 3, and an exhaustive numerical study of its efficiency is proposed in section 4. In the 
last section, this technique is applied to the "on-the-fly" correction of hyperreduced damage 
models. 

2 Problem statement and classical model order reduction for 
nonlinear problems 

2.1 Problem statement 

2.1.1 Nonlinear structural Problem 

Let us consider a structure occupying a continuous domain with boundary dfl. This struc- 
ture is subjected to prescribed displacements U D on its boundary dft u , over time interval 
[0, T]. Let u be the unknown displacement field, it belongs to the space U of kinematically 
admissible fields: 



Let IAq be the associated vector space. Prescribed tractions -F D are applied on the com- 
plementary boundary dVLf = dQ\dQ u - Under the assumptions of quasi-static evolution of 
the structure and small perturbations, the weak form of the equilibrium and the constitutive 
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relation read, at any time t G [0, T]: 

Vu* G Wo, find uGU such that: 

Jan, 




where a is the Cauchy stress tensor and e(u) is the symmetric part of the gradient of displace- 
ments. The constitutive relation between a_ and e(u) is nonlinear and described using internal 
variables (for instance damage or plasticity). It is assumed to be local and rate- independent. 



2.1.2 Space discretization 

We perform a standard finite element approximation of the space of unknown displacement 
field IA (and a similar approximation of the space of test functions Uq): 

( n n \ 

U h (n) = I u(x) | u(x) = Ni(x) Ui \ (3) 

where n n is the number of nodes, x is the position vector, (iV^jgp nn j are the finite element 
shape functions, and («Ji e [i „ n j the nodal values of the displacement field. 

The introduction of the finite element approximation into equations ([!]) and (|2j) leads to 
the following nonlinear vectorial equation at any time t £ [0, T]: 

Fmt((U| T ) ) +E Ext = (4) 

where U G M n " (n u is the number of nodal unknowns) is the vector of nodal displacement 
unknowns, F Int G M™" and F Ext G M n " are the internal forces resulting from the discretization 
of the internal virtual work (left-hand side of ^ ) and the external forces resulting from the 
discretization of the external virtual work (right-hand side of Q), respectively. 

2.1.3 Time discretization 

The nonlinear solution strategy used here is a classical time discretization scheme for quasi- 
static and rate-independent problems. This procedure consists in finding a set of consecutive 
solutions at times (tn) ng [ nt j (see [26] ) . Hence, the constitutive law ([2]) is discrete in time. 
This time discretization scheme yields the following vectorial system of n u nonlinear equations 
at time t n+ \: 

E '»'(( II l'™) me ,o„ +1 ,) +E - = !! (5) 

For the sake of clarity, the dependency of the internal forces with the history of the unknown 
fields will not be written explicitly. 
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2.1.4 Nonlinear solution strategy 



At each time i n +i the nonlinear system of equations is solved by a Newton-Raphson algorithm. 
At any iteration (i + 1) of the algorithm, the linearisation of equation §5§ around U* leads to 
the following prediction stage: 

K^<5U i+1 = -E l (6) 

where we have introduced the residual vector R = F Ext + F Int (U), and the correction of 
the increment of displacement <5U l+1 = U* +1 — U*. This linear prediction is followed by a 
correction stage: 



.1+1 _TT i XT' — ^Elnt(U) 



R l+1 = F Ext + F Int (u J+i ) k; 



=t du 



(7) 

u=u l+1 



2.2 Model order reduction 

The set of equations ([5]) can be very large if multiscale problems are simulated (see for instance 
[51 IS]). Various strategies can be used to reduce the size of this system without loosing 
accuracy. We apply here the classical model order reduction by projection to the solution of 
problem 



2.2.1 Reduction by projection 

The solution vector is searched for in a space of small dimension (several orders smaller than 
the number of finite element degrees of freedom). Let us call C the matrix whose columns 
form a basis of this space: 

Q = ( C 1 C 2 ... C n c ) (8) 

where uq is the dimension of the reduced space, and (C fc ) fcg jM nc j £ (M ra ") nc are the chosen 
reduced basis vectors. Applied to the reduction of problem (p™, the solution field is searched 
for under the form: 

U = U| in +£« (9) 

The residual of equation ^ is constrained to be orthogonal to a space of small dimension, 
which can, in theory, be different from the one spanned by (C fc )/ Cg | 1 „ c j. Let us limit our 
study to the Galerkin procedure, which writes the following constraint: 

g T R = (10) 

Hence, the reduced form of problem ([5]) reads: 

Q T (F Ext + F Int (U |tn + Q a)) = (11) 



2.2.2 Solution procedure for the reduced nonlinear problem 



At each time step of the time discretization scheme, the reduced problem (11) is solved by a 
Newton-Raphson algorithm. The (i + l) th prediction stage is performed using the following 
set of hq linearised equation: 



R 



(12) 
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where the residual of the reduced nonlinear problem (11) is defined by R fi = C T F Ext + 
C T F Int (Uu |Cq) and the increment da 



i+1 



K i+i 

=T,R 



a 1 . The correction stage reads: 
nT 3F Int (U |tn +ga) 



Sen 



"(13) 



There is a simple link between the previous solution procedure for the reduced problem 
(11) and the systematic reduction of the prediction stages performed on the full nonlinear 
problem ^ since can be expanded as 



K 



d¥ lnt (U ]tn + Q<x) 



T,R 



da 



, T 5F Int (U) 



u=u. 



-Ccci da - 



This result can be injected in equation (12), which yields: 

C T + KIC Sa i+1 ) = 



(14) 



(15) 



This last equation is a simple reduction by projection of the (i + l) th linear prediction of 
the Newton- Raphson scheme ^ used to solve the initial problem ([5]). In other words, the 
prediction stage of any Newton iteration performed on the reduced problem yields exactly the 
same solution as a reduction of the linear prediction stage used to solve the initial problem. 



2.2.3 POD reduced approximation space 

The snapshot-POD is a projection-based model order reduction technique which requires the 
knowledge of a representative family of solutions to the global problem. This set of vectors 
is called (S^/ugp „ s j and forms an operator § = ( S 1 S 2 ... S ns ). The aim is to find 
an orthonormal basis (C fc )/% g |i nc j, of dimension uq smaller than no, such that the distance 
between spaces Im(C) and Im(S) is minimum (in the sense of the Frobenius norm). 

This problem is classically solved by computing the singular value decomposition of S: 

i = M!X T = 5Zs fc u fc v fcT (16) 

k=l 

where H = ( U 1 U 2 ... U"" ) and Y = ( Y 1 Y 2 ••• Y" s ) are two orthonormal 
matrices of respective sizes n u and ns and the upper block of S is a diagonal matrix of 
positive entries (S fc ) fcg | 1 ng j (singular values), ordered decreasingly (the lower block is a null 
matrix). One can show that the optimum reduced snapshot basis (C fc )fcg[i nc j is such that: 

Q = ( u 1 u 2 ... U nc ) (17) 



2.3 POD for the reduction of damage evolution problems 
2.3.1 Damageable lattice structure 

We focus on a very simple damageable lattice structure, made of bars under traction or 
compression (figure (U])). Each of the elementary bars occupies the domain such that Q, = 
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Ube[in 6 ] ^» wnere n b is the number of bars. The direction of the bar is n kl = PkPi/\\PkPi\\ 
where Pfc and Pi are the two nodes connected to Bar b. The displacement field in Q b is 
searched for under the form: 

u = u{y)n kl (18) 

where y G [0 ||PfcP;||]. We suppose that no volume force is applied on the structure, that 
the material of each elementary bar is isotropic and homogeneous (constant material stiffness 
along the direction of the bar), and that the section of the bars is constant and equal to 
(Sb)be{in b }- Hence, the unknown displacement field is linear along the direction of the bar. 

A very basic constitutive law, based on classical damage mechanics [I], is used to described 
the behaviour of the lattice structure. The lineic strain energy of each bar reads: 

e d = l -E(l - d)S b e 2 (19) 

where e(u) = u t y(jj) is the strain measured in the direction of the bar, and d is a damage 
variable which ranges from (undamaged material) to 1 (completely damaged material) and 
is constant within each elementary bar. Two state equations can be derived from the strain 
energy. The first one links the normal force N = SbnFj.a.n^i to the strain, while the second 
one links a thermodynamic force Y to the strain: 

N= d -^ = E{l-d)S b e Y = -^= 1 -ES b e 2 (20) 

An evolution law is defined to link the damage variable to the thermodynamic force: 



d = min 




The results presented in this paper are obtained with the following set of data. The 
sections and Young's modulus are unitary. The material parameters are f) = 0.5 and a = y/2. 
Each vertical or horizontal bar has a unitary length. 

The finite element discretization of the lattice problem is done by considering that each 
bar is a linear finite element. Because of the previously made assumptions on the loading and 
geometry, the finite element solution yields the "exact" solution to the problem. 



2.3.2 Newton/arc-length solution scheme 

The family of problems that we want to solve is described in figure ([!]) . The load is applied on 
the top surface of the structure, in the y direction. The damage state represented here (light 
gray corresponds to d = 0, while darker bars are damaged, and ruined bars are removed) is 
obtained in 60 time increments, each of which is converged to a very low level of error. 

This problem being unstable, a local arc-length procedure is combined to the Newton 
algorithm. This classical \27\ |2"8] procedure is usually derived for a particular case of the 
modified Newton algorithm, namely the updated secant Newton. In our work, we need a 
certain versatility to choose the linearised operator. Therefore, the local arc- length procedure 
is detailed here for tangent and modified Newton algorithms, in a general manner. 
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Figure 1: Reference problem. The structure is a lattice made of damageable bars in traction 
or compression. The problem becomes unstable as damage localizes and propagates. 



At time t n +i, this procedure simply consists in relaxing the norm A of the prescribed 
loading, while adding a local constraint on the increments of displacement: 



-Int 



Ur, 



me[0 n+1] 



+ AF Ext = 



U 



|*T, 



m£[0 n+1] 



AU = Al 



(22) 



where F Ext is normalized, and AU = Uu_ t — Uu_ ti[ . c is a formal boolean line operator 
which extracts the maximum difference of displacement increment undertaken by an element 
(whose stiffness is strictly positive) over the current time step: 



iS[0 n+1] 



AU = max 

VAf e (lb, d\ M ,t=t n < 1 



(AU|j 



AU| Pfe ) .n kl 



(23) 



where AU|p fc and AU|p ; are the local increments of displacements over the time interval 
[t n ,t n+ i], respectively at corners Pk and Pi of the bar. 

These two equations are linearised around point (AMJ 4 ), which yields, at step (i + 1) of 
the Newton algorithm, the following linear prediction: 



5X i+1 -- 
<5U i+1 



-R. 



cont 

-d^r 1 (R^ q +^ +i F Ext ) 



(24) 



with the definition of the residuals of system ( [22] ), R cq = F Int + AF Ext , -R C ont = AZ — c AU and 
the increment 5X l+1 = A* +1 — A\ During the correction stage, the tangent stiffness operator 
and residuals are updated similarly as described for the classical Newton algorithm, equation 
([7]), while the new extraction operator is corrected as follows: 



U 



t+l 



» (Hit, 



m£[0 n\ 



(25) 



As suggested in [281 IE] , we use a quasi-Newton approximation of the previously described 
tangent Newton-Raphson algorithm. The tangent operator at iteration K^, is replaced by the 
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stiffness matrix K l obtained by setting the internal variables to previously obtained values 
(i.e.: last converged solution, or last Newton iteration). This procedure does not yield the 
quadratic convergence rate of the genuine Newton-Raphson algorithm. Yet, it is widely ap- 
plied in damage simulations for its simplicity (no need to compute a tangent operator) and 
robustness. 

The same arc-length procedure is of course applied when reducing the nonlinear problem 
by projection on a reduced basis basis. The expressions of the tangent stiffness and residual are 



similar to the one given in section (2.2.2) (the residual is slightly modified by the introduction 



of the loading parameter A). The linear prediction performed on the reduced problem reads: 



5X t+1 = 
Sa i+1 



"-^COnt _ ^* GCKrpp) 



R 



■Ext J 



(26) 



(K^ R )- 1 (R J R + a m g T F Ex t) 



2.3.3 Model order reduction 



Reduced basis 




Figure 2: Computation of the snapshot made of seven solutions obtained at the very beginning 
of the damage process (left) and reduced basis obtained by extracting the singular vectors 
associated with the three highest singular values of the resulting snapshot operator (right) 

In order to reduce the family of problems described in the previous sections by a POD- 
based projection method, we compute a snapshot of 7 solutions (loads applied on the top 
surface of the lattice, at positions x G {2 5 8 10 12 15 18}), as illustrated in figure The 
successive loads have a very low value, so that the structure is virtually undamaged. This 
snapshot is reduced to three basis vectors by the computation of a truncated SVD. 

Within the family considered to construct the snapshot, the particular problem which 
will be simulated in the following is represented in figure ([3]). The vertical load is applied at 
points x G {7 8 9}, and the time of the analysis is discretized in 50 time steps. The reference 
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solution is converged to a very low value of the norm of the residual of the Newton/arc-length 
algorithm (norm of the residual scaled by the norm of the right-hand side term equal to 10~ 6 ). 




Fi gure 3: Superposed reference solution and solution obtained with the C-POD, for i^New — 
4.10" 1 (top), additional reduced basis vectors obtained as a correction of the reduced basis 
(bottom). 



Given the normalized current solution U 



Ui 



and normalized reference solution 



U 



refi, 



ii tt re |f M , the solution error used to assess the accuracy of the current solution is: 

II— rcf ItH 2 



solution error 



u, t -u 



ref| t 



(27) 



In the following developments, a single value can be provided for this error. It is the maximum 
over time of each of these errors. 

In the very early stages of the simulation, the error is very low, which can be seen in figure 
([5j right) (the dashed line is the error defined by equation (27) as a function of time, while 
the dotted one is the value of the maximum damage in the structure). However, as damage 
increases, which eventually leads to the cracked area represented in figure ([3]), a simple linear 
combination of the modes of the reduced snapshot is obviously not sufficient to reproduce 
accurately the solution, which results in a high level of error. 



3 Error-controlled adaptive model order reduction 
3.1 Principle 

In order to handle complex changes in the topology of the structure, we propose to enrich 
the reduced basis C, "on-the-fly" , during the Newton iterations performed on the reduced 
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problem (11). Prior to the actual description of the algorithm, let us remark that, at each 



iteration of the Newton algorithm, one can compute: 



the relative norm of the residual of the reduced problem (11), i.e.: 

S — ' ||C T F 



_ .Ext II 2 

the relative norm of the residual of the initial problem ([5]), i.e.: ^ 



-Ext II 2 

Our aim is to control both these values during the solution process. More precisely, we will 
proceed to the next time increment only if the former is lower than i^New,R an d the latter 
lower than i^New, with ^New,R ^ ^Now Typically, ^New,R is set to 10~ 6 in our tests, while the 
value of i^New (which, in our tests, ranges from 10 _1 to 8.10 -1 ) will determine the accuracy 
of the successive Newton solutions, and is an important parameter of the strategy. 

The following algorithm is applied. A sufficient number of Newton iterations are performed 
on the reduced problem, so that the norm of its residual is small enough (/cR es times smaller 
than the norm of the residual of the full problem, with kji es > 10 a parameter of the method 
which can be tuned depending on the application and implementation). If the norm of the 
residual of the full problem is still higher than a desired value ^New> the following linear 
prediction is enhanced by performing a few iterations of Conjugate Gradient on the linear 
prediction of the full problem ([5]). The resulting solution is eventually added to the reduced 
snapshot basis, as detailed in the next sections. 

This process is repeated, until convergence of both the reduced problem (||Rij||2/||C T F Ext ||2 < 
^New,R), and the full problem (||R||2/1|F Ext || 2 < ^New)- 

3.2 Corrections on the linear predictions by an iterative solver 

When required by the previously described algorithm, the enhanced linear prediction is per- 
formed by initializing a projected conjugate gradient and performing a few iterations on the 
linearised problem ([6]). The stopping criterion on the normalized residual is denoted by z^cg 
and set to a high value, typically of the same order than ^ ew . The reasons for choosing this 
iterative algorithm are the following: 

• The construction of a Krylov basis requires only matrix/product operations. There- 
fore, no matrix assembly is required, which makes it particularly suitable for parallel 
computing, and ensures the versatility of our method. In the case of serial computing, 
no matrix factorization is needed, which permits to save (i) memory, as we keep the 
sparsity of the operator (ii) time if only a coarse solution is needed, which is the case in 
the following developments. 

• The projection framework is an ideal tool for the particular issue of enriching the reduced 
basis generated by the POD. Indeed, the initialization of the projected algorithm is the 



linear prediction of the reduced problem, equation (12). The complementary part of 
the solution is searched for is a space orthogonal to the reduced basis, which extends 
immediately the solution space and reduces the number of iterations required to reach 
z^cg (left preconditioning). To summarize these ideas, each iteration of the CG is a 



correction of the solution to the reduced linear prediction (12) such that this correction 
is orthogonal to the initialization. 
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3.2.1 Projected Krylov solver 

The classical projected (or augmented) conjugate gradient [29] is applied to the approximate 
solution of ^ (the linearised operator is assumed symmetric, positive and definite), which is 
recalled here: 

K T <5U = F (28) 

Where the i superscripts indicating the iteration number have been dropped for the sake 
of clarity, and F = — R. The chosen augmentation space is the previously defined reduced 
snapshot space Im(C). 

The fundamental idea of the augmented Krylov algorithms is to split the search space into 
two supplementary spaces: 

M n " = Im(g) Im(£)^ T (29) 

J_k designing the K T — orthogonality. The unknown solution vector is decomposed in this 
space under the form: 

f <5U C = CSot G Im(C) 

ffi=«Jc + fflK where | ^ = (30) 

The K — orthogonality is ensured by introducing a projector P such that: 



<5U K = P<5U K 

£ T K T P = o 



(31) 



Hence, the decomposition of the correction of the displacement increment reads: 

SV = QSa + P = SlJ K (32) 

with the classical projector: 

p = i d - g(g T K T g) _1 g T K T (33) 

This separation of the search space into two subspaces Im(C) and Im(P) in direct sum 
leads to the following uncoupled equations: 

(£ T ^ T Q)da = £ T F <5U C = Q(Q T K T Q)~ 1 Q T F (34) 

It i) = F - K T ^U c (35) 



Note that R c = F - K T ^U c = E T F and K T P = P T K T P. Equation @ is a coarse 
initialization of the projected conjugate gradient. As stated previously, it is equivalent to 
the linear prediction of the reduced problem (equation (12)). Equation (35) is the linear 
prediction of the full problem projected on Im(C)~ L =T . This system is symmetric and can be 
solved by a preconditioned conjugate gradient. Hence we solve: 

M~ 1 It E) ^-K = M ^ T 5=0 ( 36 ) 

where M is a left preconditioner (symmetric, definite and positive). In our test cases, M is 
a diagonal matrix whose entries are the elements of the diagonal of K T . Algorithm |lj describes 
the augmented preconditioned conjugate gradient. 
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Algorithm 1: Augmented preconditioned conjugate gradient 



Compute K T g, (Q^Q)- 1 ; (P = I d - g (g T K T e) 1 £ T K T ); 

<5U C = g(g T K T g)" 1 £ T E and £U Ko = 0; 

Ro = R-I T ^Uc; 

z = gM" 1 B c , w = z ; 

for j = 1, . . . , n do 

= (&,-_!, W^O/^W^j, W H ) 

= Ry-i - aj-i| T W rl 
Zj- = PM 1 R 7+1 

end 



3.2.2 Adaptation of the reduced basis 

The augmentation space, also used as reduced space for the projection of the linear prediction 
performed on the reduced problem, is chosen as: 

£ = ( £p D =NewT ScG ) ( 37 ) 

where: 

• ^pod * s ^^ia! reduced snapshot (classical POD). It is not changed, unless the size 
of the adaptive reduced basis becomes too large. This issue is discussed in the following. 

• C NcwT is a set of converged solution vectors obtained during previous time steps, when 
at least one enhanced linear prediction has been performed (otherwise the solution 
obtained is a simple linear combination of the reduced basis vectors). At the end of 
such a time step t n , the solution \J_ t=tn is orthonormalized with respect to C poD and 
C N wT by a Gram-Schmidt procedure (using the classical Euclidian dot product). This 

orthonormalized vector is denoted by Ht=t„ and used to adapt C: 

^NewT| i=tn+1 = ( ^NewT|i= trl ) ( 38 ) 

• C„„ is a set of solution vectors obtained during the Newton process at the current 



time step, let i + 1 be the current Newton iteration, the linear prediction requiring an 
enrichment by the augmented Conjugate Gradient. If <5U is the solution obtained at 
the end of this projected iterative linear solver, it is added to the reduced basis: 



£cG i+1 = ( £cg ~ ) (39) 

where S\J q is the initialization of the projected iterative solver, given by solving the first 



row of system (34). The linearised reduced operator in equation (12) is recomputed 
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i+1 

iui in^va , vv iiiuii y iuiuo u±.lc J-i.^ vv oui u uiun v v. v.. i v^i 

Note that ,by construction, (<5U* +1 — ^U c ) is K 4 — orthogonal to the previous reduced 



and the linear prediction (|12|) is performed, which yields the new solution vector U 



basis C\ 

When the convergence of the current time step is achieved, C„„ is emptied. 

The size of the reduced basis needs to be controlled. Indeed, the size of C. T _ increases 
during the simulation. Several techniques can be used when the number of additional vectors 
exceeds n c ,NcwT,Max: 

• Keep only the last vectors of C NewT . 

• Keep a few linear combination of the vectors of C^ T m , these linear combination being 
extracted by a SVD computation. 

• Consider the following operator as a Snapshot: 

l=(£ PODMn §NewT |t ^ ) (40) 

where C NewT is the set of all the previsous converged solutions, and choose as the new 

reduced basis the vectors associated to the highest singular values of S, in decreasing 
order. This new basis is split as follows: 

=|t=t„+i = ( =POD|t=t„ +1 =NewT|i=t„ +1 ) ( 41 ) 

where, CL„„. and C_„_. have the same number of columns, and C„ „, 

—POD\t=t n+ i —POD\t=t n —NewT\t=t n+ i 

have a number of columns lower than the one of C _ (divided by two in the 
following examples) . C nnn usually contains the information of the initial reduced 

— FUD\t=t n +i 

snapshot (strongly uncoupled and associated with singular values greater than 1). 

The latter technique is the one used in our examples. It has the advantage to yield a basis 
whose vectors have a very weak energetic coupling, hence granting the method stability (we 
have observed that the two former ideas, despite the orthogonality, between the vectors of 
the reduced basis, easily lead to the stagnation of the Newton process when instabilities and 
multiple solutions have to be handled). 

4 Behaviour of the adaptative strategy 

We have introduced several parameters in the corrective model order reduction strategy (de- 
noted by C-POD). The aim of this section is to provide an understanding of the sensitivity 
of the solver to a variation of these parameters, and when possible, to provide optimal or 
practical values. This study is not required to understand the method, and the reader might 



want to refer directly to the results given in section 4.3 
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Algorithm 2: Corrective model order reduction strategy (the control of the size of the 
reduced basis is not taken described) 

Load the reduced basis C ; 

Initialise the solution of the reduced problem ; 

a° = 0; 
for n = \..n% (time loop) do 

Construct the linearized reduced operator 



K i+i 

=T,R 



,T agmtCHn^l+gg) 
da 



Compute the initial residuals 

R° = Fg xt + F Int (U n _ 1+ ga ) 

Save the reduced basis 



C° <- C while 



IIR 



R\\2 



if 



IIR 



< k R 



IIRI 



IIRI 



2 



-Bret! I 2 



> ^iVeio (Newton loop) do 



2 



then 



1^ E-Ext\\2 \\£-Ext\\2 

Assemble the global stiffness 

9F Int (U) 



dU 



Save the reduced basis 

Compute the reduced basis correction by the projected conjugate gradient 

5U~-Kt 1 R* where SU = Ca + SXJ K 
Update the reduced basis 

= \\m K h 

Project the previous solution in the new reduced space 

a'^g)- 1 (g T g Tmp « 5 

end 

Solve the reduced linearized system 

<k* m = -(ly- 1 ^ 

Update the solution 



a 



i+l 



da 



i+l 



Compute the residuals 

R*+ 1 = Fg xt + F Int (U n _ 1 +ga*+ 1 ) 

r ! r +1 = g T R J+1 

i +- i + 1 
end 

Update the reduced basis if a correction has been performed ; 

Q <- ( Q° U ) ; 
Initialise the solution of the reduced problem for the next time step ; 



a" 



end 
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4.1 Error criteria 



In a first set of numerical experiments, we investigate the influence of the precision of the 
conjugate gradient iterations on the efficiency and accuracy of the method. We use the 
following set of fixed criteria: k^ es = 10 3 , nc,NewT,Max = 3. The relative norm of the residual 
of the full problem is successively set to (a) i^New = 4.10 -1 and (b) z^Ncw = 5.10~ 2 . The results 
obtained when decreasing threshold values uqg are reported in table ([!]), Experiment A. The 
following tendencies can be observed: 

• If ucg > ^New (i-e.: the conjugate gradients are solved finely compared to the actual 
maximum level of error required for the convergence of the adaptive strategy), the 
accuracy of the algorithm (i.e: the maximum value over time of the error in solution) 
increases with decreasing values of z^cg- 

• Unless z^cg > ^New (i-e.: if the conjugate gradients are solved very loosely), the value of 
z^New has no influence on the accuracy of the algorithm. 

• The number of enhanced predictions performed to reach a given accuracy is relatively 
independent on the accuracy of the conjugate gradients resolutions, hence independent 
on v C G- 

• The number of cumulated iterations of the conjugate gradients increases when the accu- 
racy of the conjugate gradients resolutions increases (i.e.: i^cg decreases), which results 
in the global computation costs to increase. 

As we want to ensure a certain consistency in the accuracy of the algorithm with respect to 
the reduced model error criterion fN ew , we set z^cg = ^New Note here that the corrective 
algorithm converges even if the convergence threshold on the conjugate gradients, z^cG; is set 
to a loose value. 

In a second numerical experiment, we try to explain the role of kn es (parameter which 
sets the maximum value of the reduced residual norm to reach before allowing a correction 
to be performed) on the behaviour of the method. The following set of parameters is used, 
^New = 1.10 , z^cg = ^New, ^c,NewT,Max = 3 and the results in terms of efficiency and 
accuracy with different values of k^ es are reported in table Q, Experiment B. The following 
remarks can be made: 

• The accuracy of the solution is not influenced by A;R es - 

• The global number of Newton iterations increases when k^ es increases, thus increasing 
the global computation time. 

This behaviour can be explained by a qualitative analysis of the influence of kn es on the 
solution strategy. During a time step where a correction of the reduced basis is required, using 
low values of kn es results in the corrective linear predictions to be performed at the beginning 
of the Newton process. Conversely, when kn es is set to a high value, a first nonlinear prediction 
is performed with the previous reduced model, and the correction is done at convergence of 
this process. As suggested by the results presented here, the former strategy is more efficient, 
as a relevant (i.e.: corrected) reduced model is constructed during the very first iterations of 
the Newton process. 
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However, we will see, in section [5| that this conclusion can be jeopardized by the actual 
implementation of the strategy. 



EXPERIMENT A (a) 




CG 


CG 


Newton 


CPU 


Solution 




corrections 


iterations 


iterations 


time 


error 


8.10- 1 


6 


82 


257 


20 


1.75 % 


3.KT 1 


8 


180 


291 


23 


1.93 % 


1.10 -1 


8 


224 


290 


23 


1.93 % 


3.10~ 2 


8 


256 


291 


23 


1.93 % 


1.10~' 2 


8 


294 


291 


23 


1.93 % 


1.10^ 3 


8 


318 


291 


23 


1.93 % 


EXPERIMENT A (b) 




CG 


CG 


Newton 


CPU 


Solution 




corrections 


iterations 


iterations 


time 


error 


8.10" 1 


42 


545 


485 


38 


0.449 % 


3.10" 1 


41 


858 


496 


41 


0.382 % 


l.KT 1 


41 


1204 


498 


44 


0.380 % 


3.10- 2 


41 


1372 


498 


44 


0.380 % 


1.10- 2 


41 


1315 


498 


45 


0.380 % 


1.10" 3 


41 


1531 


498 


45 


0.380 % 


EXPERIMENT B 


^Res 


CG 


CG 


Newton 


CPU 


Solution 




corrections 


iterations 


iterations 


time 


error 


10° 


31 


827 


384 


36 


0.894 % 


10 2 


32 


895 


410 


36 


0.887 % 


10 4 


32 


896 


453 


39 


0.887 % 


10 6 


32 


896 


500 


42 


0.887 % 



Table 1: Influence of the parameters of the C-POD on the accuracy and efficiency of the 
corrective model order reduction 

4.2 Number of additional vectors in the reduced basis 

A new set of numerical experiments are performed in order to study the relevancy of the 
adaptive reduced model created during the successive time steps. The following set of param- 
eters fixed parameters is chosen: z^New = 10 _1 , vqg = ^New, &Res = 10 3 and the value of the 
maximum number of solutions added to the reduced model, nc,NewT,Max is incremented. The 
results obtained are given in table Q. The following general tendencies appear: 

• The number of correction steps required to obtain a given accuracy drastically decreases 
as the size of the reduced model (number of additional reduced basis vectors) increases. 

• The global CPU time rapidly decreases when using an increasing, but small, num- 
ber of additional reduced basis vectors. Though, when too high a value is given for 
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^C,NewT,Max) the costs of the assembly steps becomes dominant, and the numerical effi- 
ciency decreases. 



"■C,NewT,Max 


CG 
corrections 


CG 
iterations 


Newton 
iterations 


CPU 
time 


Solution 
error 





74 


2696 


560 


60 


0.780 % 


1 


45 


1367 


435 


40 


0.773 % 


3 


32 


896 


429 


38 


0.887 % 


10 


14 


370 


409 


36 


1.02 % 


20 


14 


360 


418 


41 


1.02 % 



Table 2: Influence of nc,NewT,Max on the efficiency of the corrective algorithm model order 
reduction 

4.3 Results 

We can finally evaluate the suitability of the adaptative algoritm to the purpose for which it 
has been designed, namely the automatic correction of a reduced model. 

The previously studied parameters are fixed to the following values: vcg = ^Newi &Res = 
10 3 and nc,NewT,Max = 3. Several simulations are performed for different values of the control 
parameter z^New (interpreted as the error criterion on the reduced model). The results are 
reported in table ([3]). For comparison, the results obtained with the basic POD strategy and 
for a classical Newton algorithms on the complete system are also given. 





CG 


CG 


Newton 


CPU 


Solution 




corrections 


iterations 


iterations 


time 


error 


POD 






266 


17 


8.27 % 


^Ncw = 8.10~ ] 


4 


46 


274 


21 


4.11 % 


^Ncw = 3.10 _] 


12 


262 


334 


26 


1.61 % 


I^New = 1.10 -1 


32 


896 


429 


38 


0.887 % 


I^New = 3.10" 2 


54 


1601 


556 


52 


0.263 % 


Z/N ew = 1.10" 2 


73 


2107 


705 


67 


0.157% 


Z/New = 1.10" 3 


132 


4351 


1539 


144 


0.00424 % 


Full Newton 






824 


434 





Table 3: Efficiency and accuracy of the corrective model order reduction 

The first global observation is that the numbers of corrections and cumulated Newton iter- 
ations increases as i^New decreases. As a consequence, the CPU time increases with decreasing 
values of the control parameter. More importantly, the accuracy of the solution is controlled 
by ^New Figure ([4j left) shows that the relation between the solution error and the thresh- 
old value is almost linear (log- log scaling). Similarly, figure ([4j right) demonstrates that we 
obtain a linear relation between the accuracy of the solution and the CPU time required for 
the computation. Hence, the adaptive algorithm acts as a continuous corrective link between 
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Figure 4: Solution errors obtained with decreasing values of ^cg> for i^ ew = 4.10 -1 (left) and 
u CG = 2.10- 1 (right). 

the basic POD strategy on the one hand, and a fine solution scheme for the full nonlinear 
problem on the other hand. 

The results detailed below help to understand how the successive corrections are performed 
during the simulation. Figure ([5j left) presents the evolution of the solution error with time, 
for the different values of i^New used in this study. A "zoom" in these results is plotted in figure 
([5j right), where we focus on the curve obtained by the application of the basic POD, and the 
one obtained by the adaptive algorithm, with i^New = 3.10 -1 . The evolution of the maximum 
damage in the structure (linear because of the local arc-length control) is also reported. Each 
vertical bar point out the time steps where at least one correction of the reduced basis has 
been performed. One can see that a correction is necessary at the very beginning of the 
simulation, to take into account the specificity of the loading. The corrected reduced model 
constructed after the time step is sufficiently relevant to allow the simulation of the initiation 
of the damage in the structure, at the precision required by the predefined value of z/N e w (a 
single additional correction is performed at Time step 10). However, when the maximum 
damage exceeds 1 (Time step 20), the rate of corrective steps increases as the fast topological 
changes must be taken into account in the reduced model. 

For illustration, and in the particular case detailed in the previous paragraph, the solution 
obtained at the end of the simulation is superimposed to the reference solution in figure Q. 
The correction of the reduced basis is also plotted and shows an important deformation in 
the region where the damage localises. Elsewhere, the value of this vector is almost null, 
which means that, far away from the "crack", the displacement field is correctly represented 
by the initial reduced basis. In other terms, the correction is local, which obviously leads to 
the conclusion that multigrid or domain decomposition solvers could help to reduce the cost 
of the corrections. This issue will be discussed later on. 
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Time step Time step 

Figure 5: Error in solution for different values of the residual error of the initial problem 
^New (left). Darker curves correspond to lower values of i^j ew . The dashed curve is obtained 
by applying the basic POD and increases significantly with time. When using the corrective 
POD (C-POD), the error decreases in a monotonic manner with decreasing values of i^N e w On 
the right, the dotted line is the maximum damage over the structure. The dashed line is still 
obtained by applying the basic POD, is reasonably low in the initiation phase (no significant 
topological changes), and increases when damage propagates. The bar corresponds to the 
time steps being corrected when applying the C-POD (z^New = 3.1CT 1 ). The frequency of 
correction increases when strong topological changes occur. The dark line is the resulting 
error in solution, and is kept to a low level. 

4.4 Discussion 

We have shown in the previous section that the error-controlled POD-Krylov algorithm is 
a bridge between a full Newton/Krylov solver and the basic POD method. This algorithm 
can obviously be used to correct the reduced model of various engineering problems. Yet, in 
the case of mechanical models involving internal variables, the computation time required to 
assemble the reduced stiffness operator (evaluation of internal forces) is far too high to yield a 
computationally efficient strategy. Indeed, although the number of unknown is considerably 
reduced, the evaluation of the internal variables and internal forces are done in each element 
of the structure. In [23], it is proposes to use the hyperreduction to reduce the cost of 
these operations. In the next section, we propose to apply the adaptive algorithm to the 
hyperreduction of damage models. 

5 Krylov corrections for the hyperreduction of damage models 
5.1 Principle of the hyperreduction strategy 

The principle of the hyperreduction, which we will denote by H-POD, was first described in 
[22J, and applied to mechanical problems involving internal variables in [23]. It drastically 
reduces the cost of evaluating the internal forces required to compute the successive linearised 
operators and residuals in the Newton algorithm for a problem reduced by the POD. The main 
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idea is to consider only a few of the equations of the initial problem to find the coefficients 
associated to each of the reduced basis vectors. If the controlled equations and the reduced 
basis vectors are correctly chosen, the norm of the residual due to the equations which are not 
controlled will be small. Note that we do not aim at questioning the validity of this strategy, 
but to use it as an efficient solver for damage models, and test the efficiency of the adaptive 
reduced basis in such a framework 



5.2 Petrov-Galerkin formulation of the nonlinear system of equations 

The starting point of the strategy is the full system ([5]), in which the approximation of the 
displacement as a linear combination of reduced basis vectors has been introduced: 

F Int (U| tn +ga)+F Ext = (42) 

This system of equations is obviously over constrained. A solution to this system is obtained 
in the classical model order reduction by using the Galerkin orthogonality condition with 
respect to the reduced basis. But one can choose any orthogonality conditions leading to a 
well-posed problem. A "natural" choice can be to satisfy only a very few number of equations 



rih of (42) (called controlled equations), with > uq 



I T (F Int [U ltn + Qql) + F Ext ) = (43) 

E is a Boolean operator of n u (number of nodal unknown) rows and columns. For each 
column of E, the only non-zero value correspond to the i th controlled equation, for i G [1 n/J, 
and is set to 1. 



In the H-POD, problem (43) is required to satisfy the Galerkin orthogonality condition 



with respect to the reduced snapshot, which finally yields the following system: 

g T H T (F Int (U, tn +Q<*\) +g T H T F Ext = (44) 



5.2.1 Nonlinear solver 



The Newton Procedure described in section 2.2.2 can be applied to the hyperreduced problem 

I, 

(45) 



(44). At iteration i of the Newton algorithm the following linear prediction is performed, 

Sa l+1 



K 



T,HR- 



R 



HR 



where Rjj R = Q T ^ T (F 

Ext + F Int (U| tn + Ca 1 )). The correction stage reads: 

Rh+r 1 = g T EE T F Ext + Q T H T F Int (U |tn + Qrf +1 ) 
5ii r F Int (U |tn +g«) 



K i+i 

=T,HR 



da 



(46) 



The resulting linearized operator K* TTT1 is square, fully populated, non-symmetric and of very 

1 ,nH 

small size x n^. It is solved by a direct solver (LU factorization). 

Note that EE T is not explicitly computed. It simply extracts the internal and external 
forces of the controlled equations used to calculate the residuals and linearised operators. The 
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integration of the internal and external forces needs only be done on a reduced integration 
domain fi/,, made of all the elements connected to a controlled node (see figure 

The same developments can be done for the derivation of the arc-length procedure, which 
is used in the examples of this section. 

5.2.2 Reduced integration domain 

The choice of the controlled equations is the keypoint of the hyperreduction strategy. In [22, 
[23] , it is shown that in order all the reduced basis vectors to be "observable" , one should choose 
to control the equations associated to the nodes whose connected elements are subjected to the 
largest strain energy from the reduced basis vectors. Hence, it is recommended to perform 
a loop on the nodes and compute the average of the strain energy due to each individual 
reduced basis vector in the connected elements. Then, for each reduced basis vector, one 
chould choose a given, and small, number of nodes with the maximum average strain energy 
observed by the connected elements. Finally, the controlled equations should correspond to 
all the degrees of freedom of these particular nodes, and the connected elements define the 
reduced integration domain. 

5.3 Modified Error-controlled algorithm 

When solving the hyperreduced system by our adaptive algorithm (the resulting strategy will 
be called corrective hyperreduction, and denoted by CH-POD), certain operations become 
cost-inefficient: 

• The calculation of the residual R. 

• The assembly of the global stiffness, required to enrich the reduced basis. 

• The conjugate gradient itself if not correctly preconditioned. 
These issues are discussed next. 

5.3.1 Calculation of the damage variables and residual of the initial problem 

In the H-POD method, the calculation of the internal variables needs only be done on the 
reduced integration domain Q^. In [23], it is proposed to interpolate these internal variables 
over 0\f2/j thanks to pre-calculated snapshot functions. In our case, this particular feature is 
not of interest, for: 

• The residual of the global system needs to be controlled. 

• The residual of the global system needs to be calculated when required to perform an 
enhanced prediction. 

Therefore, the internal variables are computed systematically when the relative norm of 
the residual of the reduced problem reaches a treshold value. The norm of the full residual can 
thus be computed, and if its value is too high, an enhanced linear prediction (e.e.: a correction 
of the reduced basis by the Krylov solver) is performed. Note that this calculation is very 
expensive compared to the cost of the Newton iterations performed on the hyperreduced 
system, and must be limited to a minimum number. 
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Let us recall that in the former algorithm given in section 3.1 an enhanced linear prediction 

IIRJU IIBII2 
could be performed if — = < fcRe S ~ It requires the computation of the full 

\\Q F Ext || 2 l|l-Extl|2 

residual at any iteration of the Newton solver, which is, in the hyperreduction strategy, 
computationally inefficient. 

Hence, for this particular application, this comparison of the two norms is bypassed, and 
an enrichment of the reduced basis will not be made unless the relative norm of the reduced 
problem is lower than k'Ne-^R (convergence criterion on the reduced nonlinear problem, set to 
a very low value) , thus limiting the number of Newton iterations at which the internal damage 
variables must be calculated in f2\0/j. An another modification permits to suppress a large 
quantity of the numerical costs. The norm of the full residual is controlled at the first time 
step, and, if no correction is needed, no control is done during the following nc : T th timestep 
( n c,T = 2 in the following test cases). The same procedure is repeated after each controlled 
time step. 



5.3.2 Global stiffnesses operators 

Obviously, the assembly of the global stiffness is the weak point of this error-controlled hy- 
perreduction strategy. Yet, as the conjugate gradient algorithm is only used to enrich the 
reduced basis, one does not need to use the current stiffness matrix. The "naive" technique 
to suppress the costs associated with this operation would then be to use a stiffness matrix 
assembled during the computation of the snapshot. Though, as the damage in the structure 
increases, the relevance of the computed additional reduced basis vectors can be questioned, 
leading to a slower convergence of the adaptive strategy, which will be illustrated later on. 

We propose an intermediate technique called the "patch assembly" (PA). Before an en- 
hanced prediction stage, the residual of the full system is known. We can reasonably as- 
sume that the part of the domain £Ipa of highest local residuals correspond to high levels of 
structural changes. Therefore, we get an approximation of the current stiffness operator by 
removing the previous elementary stiffnesses corresponding to the elements in £Ipa from the 
last stiffness approximation, and replacing them by their updated values. 

This procedure is performed in three steps: 

1. Identify the nodes corresponding to the highest values of the residual. 

2. For each of this nodes, identify the connected elements, whose union defines fipA- 

^ i 

3. In the existing approximation fo the stiffness operator K , remove the elementary ma- 
trices corresponding to the elements in £lpA- 

It = St- E 4SK T ^oid)A £ (47) 

E/n E £n PA 

where Zoid represent the state variables and material parameters used to assemble K^, 
(A )£ g [ lnB j are Boolean assembly operators, which extract the rows of the global 
stiffness corresponding to the degrees of freedom of Element E, and K T (Zqw) the 
elementary stiffness operator. 
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4. Calculate the elementary matrices corresponding to the elements in and assemble 
them in K T : 

fi T fl =I T + £ AjK Ti£ (^)A s (48) 
E/n E en PA 

The performance of this simple technique will be demonstrated later on. 

5.3.3 Preconditionning 

As the size of the initial problem increases, the preconditioning of the conjugate gradient used 
to enrich the reduced basis becomes an issue. In order to be consistent with the projection 
framework that we have introduced, we make use of the Selective Reuse of Krylov Subspace 
(SRKS), proposed in [30]. Briefly, the preconditioner itself is unchanged (diagonal matrix 
whose entries are the inverse of the diagonal elements of the stiffness operator), but additional 
reduced basis vectors are added to the augmentation space: 

fi = ( £ =SRKS ) («) 
The columns of C gRKg are made of eigenvectors of the successive linearised operators projected 

on Im(C — T ) during the computation of the snapshot. These eigenvectors are the converged 
Ritz vectors obtained when solving the linear systems by a conjugate gradient (eigenvectors 
of the complete linearised operator associated with the highest eigenvalues). 

This feature permits to rapidly access the super-convergence of the conjugate gradient. 

5.3.4 Hyperreduced arc-length strategy 



The search for the largest increment of displacement in equation (23) is only done on the 
reduced integration domain. A sufficient number of point corresponding to the highest damage 
rate is set as controlled nodes at the end of each time step. Therefore, one can expect that 
the control equations is reasonably satisfied. 

5.4 Efficiency for parametric simulations 
5.4.1 Presentation of the test cases 

The test case studied here is a square damageable lattice (the constitutive modelling is the one 
given in section 2.3.1) represented in figures ([6]). The lattice structure is made of 14,520 bars 



connected by 3721 nodes. Hence the number of degree of freedom involved in this problem 
is 3721 x 2 = 7442. A non-zero homogeneous traction force such that F G = F x x + F y y is 
applied on the right edge of the structure, while zero-displacement are prescribed on its right 
edge. We define the angle: 

= arctan ( -f J (50) 



Fx, 

The load is applied through a reinforced frame occupying domain f2p. The Young's Mod- 
ulus of the bars belonging to the frame is equal to one, while their damage parameter a is set 
to 0.1. Hard inclusions, occupying domain £l\ are located in domain fi\f2p (see their position 
and shape on figure (loT)). The Young's Modulus of the bars such that f2& £ is set to 10 
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Figure 6: Reference problem and associated solutions obtained for different material param- 
eters and direction of the prescribed load (left), where darker colours correspond to higher 
value of damage, and Young's modulus cosinusoidal distribution in these two cases (right), 
where darker colours correspond to higher values of the stiffness parameter. 

while a = 0.01. Due to these low values of a in Vt-p U fii, the damage will localize in the 
complementary domain Qq = f2\(Jlp U Oi), where a is set to 1. The young modulus of the 
bars in this domain is defined by the following law: 

E\ Qbe n c =E C (1 + e c (sin (w c ((X b - X) + (Y b -Y)))+ sin (u c {(X b - X) - (Y b -?))))) 

(51) 

where the position of the barycentre G b of a bar occupying domain is 

OG b = X b x + Y b y (52) 
and the position of barycentre G of the structure is defined by: 

OG = Xx + Yy (53) 
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The purpose of our study is to obtain the dissipated energy in the structure, for a value 
of the load that is close to the instability point (hence the arc-length control is not required) , 
as a function of the two following parameters: 

• the angle (ft, which controls the direction of the force. 

• the angular frequency ujq, which controls the spatial distribution of the Young's modulus 
of the bars belonging to 

The construction of such functions is very classical is stochastic analysis. It is often associated 
with the Monte-Carlo method to determine the random distribution of a quantity of interest 
(in our case the dissipated energy) as a function of the stochastic input parameters (here (j> and 
ojq)- If an analytical derivation of this function is not possible, the classical procedure consists 
in performing a few deterministic simulations, for different values of the input parameters, 
and interpolate the value of the quantity of interest between these points. 
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Figure 7: Estimation of the response surface. The surface relates the dissipated energy to 
two structural parameters (j> and uq over a rectangular domain. The values of this function 
at 25 nodes (crosses) are computed by performing nonlinear simulations. The other values 
are interpolated using a piecewise linear approximation. 



Let us recall that the dissipated energy of the chosen damage model is given by: 
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The following values are chosen for the constant parameters: Eq = 1, ec = 0.2. Particular 
distributions of the Young's modulus are represented on figure ([6j right). 

We will consider that the input parameters can vary within the following range: 



w c G [wc,i,wc,2] = [0.05,0.1] 

<l>e[<h.,<h] = [0,5] 



(55) 
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25 simulations are performed to estimate the response surface (see figured 



f uj C G {0.05,0.625,0.75,0.875,0.1} , . 

\ G {0,1.25,2.5,3.75,5} ^ j 

The response surface (linear interpolation) obtained when using a direct Newton solver over 
10 time increments is given in figure ([8j left). The displacement solutions, damage maps and 
Young's modulus distributions corresponding to the two limit points boxed in figure ([7]) are 
represented in figure Our goal is to limit the costs of these nearby simulations by using 
POD-based reduced order model for these 25 successive solutions, and in particular to test 
the efficiency of the CH-POD. 



Basic H-POD 




Figure 8: Response surface obtained in the reference case, and when using the H-POD. The 
error in the dissipated obtained when not correcting the reduced model is very high, due 
to the large distance between the solution which is searched for and the a priori computed 
snapshot. 

5.4.2 Parameters of the hyperreduced model 

In order to build a snapshot, we solve successively four problems circled in figure [7} Each of 
these simulations is performed in 10 time increments. The final snapshot, made of 40 vectors, 
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is reduced to 8 reduced basis vectors by a singular value decomposition (the ratio of the first 
singular value of S over the eighth is 1.9 10 -8 , which shows that the snapshot is correctly 
represented in C). 

We construct the reduced integration domain by choosing a list of controlled nodes (2 
controlled balance equations per node). 

1. The domain is divided into 100 square subdomains. For each of these subdomains, an 
arbitrarily chosen node is set as controlled point. 

2. 5 nodes are arbitrarily (first ten in the global nodal numbering) chosen for each of the 
prescribed boundary conditions (left and right edge of the lattice). 

3. For each vector of the reduced basis, 5 nodes whose connected element have the highest 
strain energy are set as controlled nodes. 

4. 20 additional nodes with the highest damage rate are added to the list, hence allowing 
to precisely integrate the constitutive law in this particular elements. 

An example of RID obtained during the CH-POD process is given in figure d9|). 
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Figure 9: Reduced integration domain (dark bars) obtained at the end the first simulation 
performed to estimate the response surface. During the Newton process, the integration of the 
internal forces need be done only on this domain, which drastically reduces the computation 
costs. 



5.4.3 Results 

The following solvers are to estimate the response surface: 

• Newton algorithm on the complete system, the linear prediction being solved by a direct 
solver, and z^New = 10~ 6 . 

• "basic" hyperreduction (H-POD) with the snapshot and controlled points defined pre- 
viously. 
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• corrective hyperreduction (CH-POD), where fNew.R = 10~ 6 (threshold value on the 
relative norm of the residual of the reduced problem), nc,NewT,Max = 3 (number of 
maximum additional vectors added in the reduced basis during the Newton iterations) 
and i^New successively set to (a) 1CT 1 ,(b) 3.10 -2 and (c) 10~ 2 (error criterion on the 
reduced model). 

For each of these solvers we report the total CPU time required to construct the surface, 
and the maximum error in the dissipated energy over the 25 simulations (see table Q). This 
latter error is defined as follows: 

Dissil Reference Dissil Current experiment t r r-j\ 

error = — (57) 

DlSSi| Reference 





Maximum error 


CPU time 




in dissipation 


00 


Full Newton (reference) 




3227 


H-POD 


113 % 


618 


CH-POD (a) 


5.06 % 


730 


CH-POD (b) 


1.98 % 


890 


CH-POD (c) 


0.82 % 


1075 



Table 4: Error and CPU time obtained when estimating the response surface with the H- 
POD and CH-POD compared to a reference response surface obtained by a classical Newton 
algorithm on the complete system of equations. 

The results given by the hyperreduction are not satisfying. The maximum error obtained 
exceeds 100%, which can be explained by observing the results presented in figure Q. One 
can see that the values of dissipated energy obtained for values of ujq which have not not been 
taken into account in the snapshot are wrongly predicted, and especially for uq = 0.05. In 
this latter case, the structure is a lot weaker (greater quantity of bars with Young's modulus 
lower than Eq, as appears in figure ([6])) than in the cases computed to create the snapshot. 
As a consequence, the dissipated energy is overestimated. 

The application of the CH-POD permits to retrieve a correctly estimated response surface 
for small additional costs. For instance, in experiment (a) (i^New = 10 _1 ), the maximum error 
drops to 5% for 18% of additional CPU time compared to the basic reduced order simulation. 



The response surface obtained is given in figure (10). 



Decreasing values of the reduced model error criterion (z^New) results in a decreasing value 



of the error, while slightly increasing the additional CPU costs (see figure (10) and table Q). 
In experiment (c), the difference between the response surface obtained and the reference one 
can hardly be notice, and the evolution of the quantity of interest with respect to the input 
parameters is correctly taken into account. Yet, this simulation is less than two times as 
costly as the one made using the H-POD. 

The CPU time, corresponding to each of the 25 simulations, in case (b), is also reported in 



figure (10). The additional costs are due to the solution which is searched for to be far away, 
in some sense, from the snapshot. As explained previously, the Young's modulus distribution 
leads to the snapshot being irrelevant, which explains why the costs of the corrections in 
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Figure 10: Response surface obtained in the reference case, and when using the CH-POD for 
various values of the reduced model error criterion, and CPU time for each of the 25 simula- 
tions (additional CPU times with respect to the ones obtained by applying hyperreduction). 
The error in the dissipated energy monotonically drops when the criterion decreases. The 
additional CPU time is concentrated in the simulations whose solution is far away from the 
snapshot. 
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the CH-POD increases. Conversely, the variation of the angle, for <p > is correctly taken 
into account in the uncorrected model order reduction, and the CH-POD generates very few 
extra costs in these cases. In the particular case where <j) = the solution is not correctly 
represented by a linear combination of the initial reduced basis vectors. Indeed, a strictly 
positive angle in the direction of the load leads to the localization of damage in one of the 
corner of the square structure (see figure [6]), while the structure is much stiffer for = 0. The 
CH-POD automatically corrects the resulting reduced model. 



5.4.4 Brief validation of the Patch-assembly 

In order the structural changes to be significant in the structure, Experiment A is performed 
with a higher value of the norm of the final loading (same value of the arc-length, but greater 
number of time increments). The evolution of the damage leads to the apparition of cracks 
in the vicinity of two of the hard inclusions. 

We use three different approximation of the global stiffness operator in the ECH-POD 
strategy. The threshold value i^New is set to 10" 1 . In the first case, the stiffness is updated at 
the beginning of each time increment, which, as described previously, leads to an inefficiency 
due to the costs of the global assembly. In the second case, the initial stiffness is used to 
perform the enhanced linear prediction. The CPU time slightly decreases but the number of 
conjugate gradient calls required to obtain a relevant reduced basis increases significantly. 




Figure 11: Solution obtained using the arc- length procedure after the limit point (left). On 
the right, the black bars are connected to the nodes corresponding to the highest values of 
the residual of the complete nonlinear system. They are mainly located in the zone of high 
topological changes. 



Finally, the Patch-assembly technique is used to yield successive approximation of the 
global stiffness matrix. The number of elements in £1 pa is 2440 (in comparison with the total 
number of bars: 14250), which gives a decrease of factor 6 of the costs of the global assembly. 
Figure ( 11 ) shows Domain £Ipa, which is indeed mainly located in the zones of high structural 
changes (the two cracked regions). The resulting CPU time is further decreased compared 
to the second experiment, as shown in table [5] Note that the number of conjugate gradient 
calls does not increase compared to the one obtain when updating systematically the stiffness 
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operator, which shows that the approximation of this matrix is satisfying. 





CG 


CG 


Newton 


CPU 




corrections 


iterations 


iterations 


time 


Updated stiffness 


23 


424 


414 


421 


Initial stiffness 


33 


500 


412 


407 


Patch-assembly 


22 


440 


447 


335 



Table 5: Efficiency of the Patch- assembly 



5.5 Discussion 

5.5.1 For further efficiency 

A few important issues can severely limit the efficiency of the CH-POD derived in this section. 

• The efficiency of the H-POD itself can be put into question for softening behaviour. 
Indeed, its numerical complexity is directly dependent on the size of the reduced inte- 
gration domain. We have observed that choosing this domain in order to ensure the 
stability of the algorithm was not a simple issue, and chosen to use a large number of 
controlled nodes. By a better understanding of this behaviour, one should be able to 
reduce the size of the reduced integration domain to a minimum, and therefore diminish 
the costs. 

• The reduced integration domain is updated at each correction of the CH-POD, for a 
new reduced basis vector needs to be "observed" . This operation is costly as it requires 
to evaluate quantities on the whole domain. This operation needs to be optimized. 

• The price of the correction is still high. As a result, if the solution which is searched 
for is far away from the initial snapshot, the method becomes computationally more 
expensive than a nonlinear direct solver on the complete system. It is thus limited, 
up-to-now, to the correction of reduced models for the resolution of nearby problems. 
In a forthcoming paper, we will propose a multilevel approach in order to perform the 
correction only in the zones of the structure undergoing the largest topological changes. 
This strategy should lead to a significant gain in terms of numerical efficiency, and 
especially in the case of localized nonlinearities, as suggested by the results provided by 
the patch-assembly procedure. 

5.5.2 Extension to finite element problems 

The proposed method has only been validated on damageable lattices so far. Slight differ- 
ences in the efficiency results might appear when dealing with 2D or 3D softening problems 
discretised by the finite element method: 

• the conditioning of lattice problems is usually better than the one obtained when using 
finite element discretisation. We have been able to use a very simple preconditioner 
for the iterative corrections of the reduced basis: deflation orthogonally to a few Ritz 
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vectors by the SRKS method, and diagonal preconditioner on the resulting deflated 
problem. Obviously, an enhanced preconditioner should be preferred when dealing with 
finite element problems, in order to keep the costs of the corrections low. 

• the connectivity in lattice structures can be very high. As a result, the reduced integra- 
tion domain associated to each of the controlled nodes in the hyperreduction method 
is quite large, especially compared to the ones obtained when using structured finite 
element meshes. Therefore, the cost reduction associated to this particular step of the 
CH-POD should increase when decreasing the connectivity of nodes. 

One of the major issues to be addressed when dealing with FE softening problems is the 
question of localization. It is well known that nonlocal damage models should be used to 
avoid spurious dependency of the solution to the mesh. The efficiency of the hyperreduction 
technique, which focuses on local patches in the structure, has not yet been investigated in 
this framework. Its extension does not seem trivial. 

5.5.3 Suitability to parallel computing 

The extension of our strategy to parallel computing is quite straightforward. The H-POD 
can be implemented in parallel as follows: the domain is divided into subdomains, assigned 
to separate processors. Solving the small problems resulting from the linearisation of the 
reduced problem is done systematically on every processor. The evaluation of the internal 
forces is a local operation performed on the reduced integration domain, distributed on the 
processors. It only requires an assembly operation, of small numerical costs as involving only 
degrees of freedom associated to the controlled nodes. 

The corrections performed in the CH-POD are very easy to parallelize as they are based 
on a conjugate gradient algorithm. The projection operation itself has been widely extended 
to parallel computing (see [9] for example). In addition, using the condensation method, 
as proposed in the Schur-complement-based domain decomposition methods would provide 
better preconditioner than the one used in this work. 

6 Conclusion 

We developed an efficient corrective tool for the adaptive model order reduction of highly 
nonlinear mechanical problems. The novelty of the approach is that it completely integrates 
the corrections inside the projection framework classically used in model order reduction. 
More precisely, the corrections are performed using an iterative conjugate gradient, which is 
enabled by the strong link between, on the first hand, POD-based model order reduction and, 
on the other hand, projected Krylov algorithms. The resulting method establishes a bridge 
model order reduction and classical Newton/Krylov solvers, ensuring its versatility. 

The robustness of the strategy has been demonstrated by correcting the reduced model 
associated with a very complex case of damage propagation. In terms of raw CPU time, our 
algorithm had been coupled to the hyperreduction method in a straightforward way, making 
it computationally efficient in the case of mechanical problems involving internal variables. 

Yet, the corrections are still expensive compared to solving the initial reduced model. 
As explained briefly in the last section of the paper, the corrections need not to be done 
everywhere in the structure, especially when localization phenomena are involved. This can 
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lead to prohibitive costs if the solution that is searched for is distant from the snapshot. 
Therefore, a possible enhancement to increase the efficiency of the adaptive model order 
reduction is to use a multilevel algorithm to perform the corrections. 
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